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In this study, we investigate the so called carbuncle phenomenon by means 
of numerical experiments and heuristic considerations. We identify two main 
sources for the carbuncle: instability of the Id shock position and low numerical 
viscosity on shear waves. We also describe how higher order stabilizes the Id 
shock position and, thus, reduces the carbuncle. 


1 Introduction 


For the simulation of shock structures in multidimensional gas flows there are essentially 
two major approaches: shock fitting and shock capturing. The idea of the first class of 
schemes is to exactly detect the shock position and split the computational domain along 
the shock line into two (or for more complex structures even more) domains, leading to an 
almost perfect reproduction of shocks in the numerical solution |TTj^ 42]. Disadvantages 
of this method include the difficulty to deal with shocks which unexpectedly evolve in the 
domain and a practical restriction to certain shock structures. In order to overcome these 
restrictions, more and more scientists started to build their simulations on shock capturing 
schemes, which are designed to work without the knowledge of the exact shock position. 
Thus, the shock is captured in a (hopefully) thin layer of grid cells. For an overview over 
shock fitting schemes and some recent developments in that area, we refer to [|6 52]. 
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A main ingredient of shock capturing schemes are so called Riemann solvers: numerical 
(first order) flux functions, which are based on an approximate solution of the Riemann 
problem at the cell face. This class was later on expanded, starting with [54 56], through 
introduction of Flux-Vector-Splitting schemes, which are based on a splitting of the physical 
flux function, but in a wider sense can be considered as Riemann solvers. More important 
for our study however, is the distinction between complete and incomplete Riemann solvers. 
While the first are designed to resolve all waves present in a Riemann problem, the latter will 
neglect some waves. Prominent examples of complete Riemann solvers are the schemes by 
Roe [50] and Osher [44]. While these solvers are preferable when complex wave structures 
as well as entropy and shear waves are expected, incomplete solvers are known to be 
robust in situations dominated by strong shocks. An example is the HLL-solver [21 ] whose 
construction is based only on the two outer waves of the Riemann problem. In the eighties 
and nineties, more and more applications were treated with above mentioned methods for 
gas dynamics. The methods were also extended to other hyperbolic conservation laws like 
shallow water or compressible magnetohydrodynamics (MHD). For a detailed discussion of 
shock capturing, we refer to the textbooks [ 16, ^ 35 ^ 61 ]. 

In the context of shock capturing, some irregularities were observed: properties of the 
discrete solutions which were by no means representations of physical phenomena. In gas 
dynamics simulations unphysical discrete shock structures and even a complete breakdown 
of the discrete shock profile could appear [47 48]. According to its form in blunt body 
problems, it was christened carbuncle phenomenon. Since the seminal paper of Quirk [48], 
an immense amount of research has been conducted on this instability problem. The origin 
of the name comes from the fact that in strongly supersonic flows against an infinite cylinder 
simulated on a body-fitted, structured mesh the middle part of the resulting bow shock 
degenerates to a carbuncle-shaped structure. It was conjectured already by Quirk [48] 
that this phenomenon is closely related to other instabilities such as the so-called odd-even- 
decoupling encountered in straight shocks aligned with the grid. Unfortunately, the failure is 
only found in schemes with high resolution of shear and entropy waves, so called complete 
Riemann solvers, which are needed to properly resolve the boundary layers and turbulent 
structures. This category includes for example the Godunov, Roe, Osher, HLLC and HLLEM 
schemes [ 13 58]. These schemes are preferable in calculations involving complex wave 
structures as well as boundary layers. 

The research on the carbuncle was twofold. On the one hand, the stability of discrete 
shock profiles was investigated in one as well as in several space dimensions. On the other 
hand, a lot of effort was put into finding cures for the failure of some schemes in numerical 
calculations. For example, many cures that were offered, are based on an indicator that 
tells the scheme when to switch to an incomplete Riemann solver. These indicators need 
information from other cell faces, making the numerical flux function non-local. It was 
found that even in one space dimension there are some instabilities of discrete shock 
profiles: slowly moving shocks produce small post-shock oscillations [|T 23,48]. But also 
in the case of a steady shock, instabilities can be found depending on the value of the 
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adiabatic coefficient y as was shown by Bultelle et al. [|7j]. However, the connection to 
two-dimensional instabilities is still not fully understood. 

The two-dimensional instabilites themselves seem to be closely related to each other. 
Chauvat et al. [1^ show through an ingenious numerical investigation that the mechanisms 


driving the odd-even-decoupling and the carbuncle are closely related. Dumbser et al. [ 12] 
present a method to test Riemann solvers for their tendency to odd-even-decoupling. Here, 
the basic idea is to discretise a steady shock in space and test the linear stability of the 
system of ODEs resulting from the Method of Lines. This allows for all tested solvers 
to predict whether they would evolve an instability or not. There is also a number of 
experimental studies of the carbuncle, especially the influence of the underlying numerical 
flux function [32 -34 59], with the goal to identify the “optimal Riemann solver”. Finally, 
Elling [14] found a connection to physical shock instabilities. 

Most these investigations have in common that they (a) intend to find a single source for 
the carbuncle, (b) do not take into account the influence of the order of the scheme—they 
usually compare different Riemann solvers in a scheme with fixed order—, and (c) do not 
distinguish between the contribution of entropy or shear waves to the carbuncle. The most 
surprising is case (b) since it is well known that in higher order schemes the carbuncle is 
much weaker than in first order; for very high orders, it is essentially absent. The purpose 
of this paper is to fill these gaps in research. We want to study the influence of the (Id) 
stability of the shock position and the 2d or 3d features such as vorticity separately. In this 
course, we also try to separate the influences of entropy and shear waves. But the main 
focus (and main novelty) of this study is that we investigate the influence of the order of 
the scheme on the stability of the (Id) shock position. We will show how increasing the 
order of the scheme, despite of lowering the numerical shear viscosity, stabilizes the Id 
shock position. 

The outline of the paper is as follows: In Section]^ we give a short representative review of 
some theoretical results. The insight gained by these results provides us with the guidelines 
for our numerical experiments. In Section]^ we give a review of the schemes we use in 
our numerical experiments. The numerical test cases are introduced in Section The 
main results are presented in Sections [^(one-dimensional issues) and [^(multi-dimensional 
issues), followed by some conclusions and directions for further research in Section [^ 


2 Short review of the theory 


There are many papers dedicated to the carbuncle phenomenon 


14 29 31 43 45 


46 48,53,63 - 65 ], however only few of them discuss the origins of the carbuncle from a 


theoretical point of view. Here, we give a short representative review of some theoretical 
results. 
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2.1 Contribution by Buiteiie et ai. 

Bultelle et al. []^ investigate steady shocks in one space dimension. A first study of these 
was done by Barth [Q] who found that for a perfect gas with y = 1.4, flux functions which 
enforce the Rankine-Hugoniot condition at discontinuities may have transition states which 
are unstable to perturbations when the preshock Mach number is greater than six. Bultelle 
et al. go even further. They prove that the Godunov scheme for strong steady shocks and 
an adiabatic coefficient 1 < y* ^ 1.62673, with y* being a root of 

/ + 3y^ - 21y^ + 17y + 8 = 0, 

can produce unstable shock profiles. They report that in practice, after a transient regime, 
the unstable leads to a stable profile with the intermediate state in a neighbouring cell or 
even one cell further. If in neighbouring grid slices normal to the shock front, the shock 
position jumps in a different direction, say +2 cells in one and —2 cells in another slice, 
this leads to an unphysical crossflow along the shock. We discuss this situation in more 
detail in Section IsTTl 


2.2 Contribution by Roe and Zaide 


while Bultelle et al [j^ discuss steady shocks in general. Roe and Zaide [63 -65 ] focus on the 
long-time behavior. Although they mainly discuss the standard Roe solver, the discussion 
relies on the fact that the solver tries to establish the exact Rankine-Hugoniot condition at 
single shocks. Thus, we can expect that their results also apply to other Riemann solvers 
with this property, e. g. Godunov, HLLE etc. Roe and Zaide investigate the behavior of 
steady discrete shocks with one intermediate point in Id. For the Euler equations (and 
also for the shallow water equations) it is impossible to find for a non-trivial steady shock 
with left and right states q; and q^, which are related via the Rankine-Hugoniot condition, 
a middle state which is related to both states via the same condition. As a result, the 
scheme enforces a middle state which is unphysical. In some situations, is not even 
constant in time: the shock, although unsteady, is trapped in a single grid cell. And even 
if is steady, the shock position in the cell is ambiguous. Depending on the conservative 
variable used to compute the shock position one gets different results. Another feature of 
these shocks is an overshoot in the momentum, which is also the mass flux, again indicating 
the deviation of the discrete solution from real world physics. 

From these results we can draw the conclusion that, at least for Riemann solvers based 
on the Rankine-Hugoniot condition, even shocks which remain in the same grid cell are 
highly sensitive to perturbations. 


2.3 Contribution by Eiiing 

Filing [ fT^ investigates the influence of the supersonic upstream region on the shock profile. 
He models the interaction of a vortex filament with a strong shock. For this purpose. 
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Figure 1: Model situation for physical carbuncle. 


he starts with a steady shock. In the region upstream of the shock, he picks the middle 
slice of the computational grid and artificially sets the velocity to zero as sketched out in 
Figure This closely corresponds to the experimental and theoretical investigations by 
Kalkhoran and Smart [24] and by Zhang et al. [66]. It turns out that the carbuncle-like 
structure which comes with the Godunov-scheme is similar to the experimental results. 
Furthermore, Filing gives numerical evidence that the structure does not depend on the 
energy inherent in the fllament. Even more so, the structure is the same for all reasonable 
resolutions of the computational domain. The numerical results obtained with Osher or 
HLLEM are almost identical to those obtained with the Godunov scheme in contrast to the 
HLLE and other incomplete Riemann solvers. We And that in real world simulations of 
a shock interacting with a vortex fllament, the viscosity, especially the shear viscosity, of 
HLLE and other incomplete Riemann solvers is too high. It outweighs the physical viscosity 
by far. This explains why Filing not only states that the carbuncle is incurable but also that 
the carbuncle should not (completely) be cured, or better: not completely be prevented. 
The challenge is to avoid unphysical carbuncles and, at the same time, allow physical 
carbuncles. In other words, one has to look for schemes which allow shear and entropy 
waves to be well resolved but still prevent the unphysical breakdown of shock waves. 


3 Review of the numerical schemes for our 
experiments 

In order to perform meaningful numerical experiments, one has to choose a suitable set of 
numerical methods. Since it is well known that the choice of the Riemann solver has a strong 
impact on the tendency of the scheme to develop a carbuncle, we have to compare Riemann 
solvers which differ from each other in some features coincide in other features. Motivated 
by the discussion in the previous section, we choose solvers with different approaches at 
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shocks: Osher as an example of a solver which abandons the Rankine-Hugoniot condition 
as a construction principle and some members of the Roe/HLL-family which in general force 
the Rankine-Hugoniot condition at single shock^ Since the carbuncle is characterized 
by strong shear flows among the Roe/HLL-type solvers, we select HLLEM, HLLEMCC, and 
HLLE in order to compare different shear viscosity mechanisms and their relation to the 
carbuncle. Eor the discussion of the influence of the order of the scheme on the carbuncle, 
we employ several standard 1st- and 2nd-order variants as described below. 


3.1 Basic finite voiume scheme 


For shallow water equations, we employ CLAWPACK [10], which is an implementation 


of the wave propagation approach. Thus, limiting for higher order is always done on 
characteristic variables. Higher order calculations are done (1) on Cartesian grids with the 


Superpower limiter [26 27 ] and (2) on non-Cartesian grids with the Albada 3 limiter 


Both are smooth limiters. While on Cartesian grids the more compressive Superpower 
limiter could be applied, on non-Cartesian grids one has to go back to a less compressive 
limiter. In Section we will discuss the reason why it is advantageous to opt for a more 
compressive limiter. The Superpower limiter is a generalization of the Power limiter by 


Serna and Marquina [55]. In contrast to the original Power limiter, the powers are not 


fixed but adapted to the local CFL-number in such a manner that the resulting limiter is 
always TVD and approaches third order behavior in smooth regions. The Albada p limiters 
are derived from the van-Albada limiter in the same way as Serna and Marquina derived 
the Power limiters from the van Leer limiter. While Power 2 is just the original van Leer 
limiter, Albada 2 is just the original van Albada limiter. Albada 3 is the most compressive 
version that still makes up a CFL-number independent TVD limiter. 

For the Euler results, we employed Euler2d, a 2d-Cartesian code developed in the Group 
of Claus-Dieter Munz at Stuttgart University. The code implements standard finite volumes. 
For higher order, direction-wise geometric limiting with minmod on primitive variables is 
used. 


3.2 Roe, HLL and their relatives 


Since the Roe scheme and the HLL-type schemes are closely related [13 28], we treat them 


as one family of schemes. For a better understanding it is convenient to start with the Roe 
scheme. 

3.2.1 Roe 


The core of the Roe scheme is a consistent local linearization [ 50 ], which was first announced 


by Roe and Baines [51 ] as a means to employ their TVD-limiters in a wave-wise manner. In 


^ We will also present a modification of HLLEM that yields a behavior at steady shocks close to the Osher 
scheme. 
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( 1 ) 


the linear case, one finds for the flux function 

nq) = Aq 

and, thus, for the flux difference 

fiQr)-fiqi) = Aiq,-qi) . (2) 

For a local linearization to ensure local conservation, it is crucial to satisfy a similar condition. 
Furthermore, the linearized system should be h 5 ^erbolic, and the system matrix should 
depend continuously on the left and right states qi and qr,.. Therefore, Roe came up with 
the following conditions for a consistent linearization with system matrix A{qi,q^): 

fiqr)-fiqi) = Mqi,qr)iqr-qi), (3) 

Mqi,qr)^Aiq) for iqi,qr)^iq,q), (4) 

A{qi, q^) is diagonalizable for all qi,qy . (5) 

A matrix A{qi,q^) that satisfies these conditions is called a Roe matrix or a consistent 
local linearization for the underl)fing system of conservation laws. If there exists a single 
state q = q{qi,qr) with 

Mqi,qr) = A{q), (6) 

then it is called a Roe mean value for qi,qy. For the Euler equations of gas dynamics a Roe 
mean value is given by 

P PlPr ’ 


U = 


V = 


W = 


H = 


^/Pl + 
/pjvi + 

V Pi V Pr 

/pjwi + 

/Pi + 
/pJHi + 


(7) 


^/Pl + 

with = u^ + v^ + in the full three-dimensional case. For 2d, we just have to omit the 
values for the third velocity component w. Similarly, we find a Roe mean value for the 
shallow water equations by 


h = -(h, +h,). 


u = 


V = 


y/lllUi + -v/h^Ur 
y/Jll + 
y/FiVi + y/vv, 

V^i + Vk 


( 8 ) 
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Together with wave-wise application of standard upwind, this results in the numerical flux 
function 

Ql) = + fiQr)) “ ^ l-4(g)| Aq , (9) 

where the absolute value is applied to the eigenvalues of the matrix. This numerical flux 
is usually referred to as the standard Roe solver and is an example for a complete flux. 
Note that the wave-wise application of schemes other than standard upwind would only 
affect the eigenvalues of the viscosity matrix. In this study, we do not use the solver in 


the form @ but with the so called Harten-H 5 mian fix [20], which slightly increases the 
numerical viscosity at sonic rarefaction waves and thus prohibits the sonic glitch. The sonic 
glitch would otherwise lead to a representation of sonic rarefaction waves as rarefaction 
shocks. It is also possible to resemble the following HLL-type schemes by simply modifying 
the eigenvalues of the viscosity matrix in 


3.2.2 HLLE 


As an example of an incomplete flux, we employ HLLE. In [21 ], Harten, Lax, and van Leer 
present and discuss a variety of numerical flux functions, the simplest and most robust of 
which is usually called HLL, a scheme with very low computational cosl|^ Their basic idea 
is to start from conservation. First, one estimates < 0 < for the bounding speeds 
of the Riemann problem given by left and right states q^, q^. If one uses conservation for 
rectangle [S;, S J x [0,1] in space and time, the mean value of the conserved quantities q in 
the intermediate states of the Riemann problem can be computed. From this, by integration 
over [S;,0] x [0,1] and [0,5^] x [0,1] and averaging, one obtains the numerical flux 

OHLLiQr, Ql) = \ fiQl)) + - Ql) ■ ( 10 ) 


We assume now that A = A{qi,qj.) is a consistent local linearization, a so-called Roe 
matrix, and thus satisfies condition @. Then ([T^ can be rewritten as 

OuLhiflr, Ql) = \(fiQr) + fiQl)) - \ “ Ql) + c i^r “ Ql) ■ (H) 


Hence, the viscosity matrix of the HLL-flux is 

Sr + S 


V = J" ' f A-2 I 


Sr Si 


Sr Si 


( 12 ) 


and has the same left and right eigenvectors as the Roe matrix A itself. The eigenvalues 
and thus the wave-wise viscosity coefficients are 


Sd + St 


A ,-2 


SpS 


R'-’L 


( 13 ) 


^Their more elaborate solvers somehow anticipate HLLC [ 57 ] 
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with Afc being the eigenvalues of A. Apparently, the choice of the bounding wave speeds is 
crucial for the numerical viscosity. In the literature, many choices of the wave speeds 5^,5^ 


are given. For an overview, the reader is referred to [ 58 ]. Here, we mainly rely on the choice 


suggested by Einfeldt [13]: If it is affordable to compute the leftmost and the rightmost 


wave speed of a consistent local linearization and and the flux function / is convex, 
then set 


Si = min{Ai, 0} , 5^ = max{A,„, 0} . 


(14) 


This ensures that for both a single discontinuity and a single rarefaction wave the estimate 
of the maximal and minimal wave speed is sharp. The resulting numerical flux is called 
HLLE. Like the standard Roe scheme, it tries to force the Rankine-Hugoniot condition at 
shocks. 

3.2.3 HLLEM 

The first scheme that exploited the relation between Roe and HLL t 5 ^e schemes is the 


HLLEM scheme for gas djmamics by Einfeldt [ 13 ]. It is an attempt to formulate the Roe 


scheme as correction to HLL. There are several advantages: The computational effort is 
reduced, the adjustment of the viscosity on the acoustic waves can be easily applied by 
choosing the bounding speeds and S^, the sonic point glitch can be avoided, and the 
failure of the standard Roe scheme for strong rarefaction waves can be healed. 

The construction is as follows: Eor the sake of simplicity, we assume for the eigenval¬ 
ues < ^2 < • • • < of the Roe matrix A that A^ < 0 < Thus, we can choose Si = A^ 
and Sn = A^. With this setting, the viscosity matrix of the Roe scheme can be written as 


Vr.o, = V„^ + ^i^RKL 

A™Ai 




(15) 

(16) 


with the anti-diffusion-matrix 


K = diag(0, 52 ,..., 5 ^_ i , 0 ) . 


(17) 


For the two-dimensional Euler equations (m = 4), we find for the so-called anti-diffusion- 
coefficients 

|u| 


= 53 = 2(1 


;)■ 


c + |u| 

The special structure K allows us to express the standard Roe flux by only using the 


eigenvectors corresponding to the entropy and shear wave. Park and Kwon [46] show 


that, independent of the choice of Si and S^, HLLEM resolves single contact waves exactly 
when adhering to the Roe mean values for the anti-diffusion-term, i. e. if we stick to ( [T6] ) 
instead of ( jls] ) as originally suggested by Einfeldt. Eor HLLEMCC, our modification of 
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HLLEM as discussed in the next section, we nevertheless employ the original Einfeldt 


setting ( [TS] ). The loss in resolution of the scheme is rather small. An advantage of ( [TS] ) is 
that it deactivates the anti-diffusion terms automatically for full upwind, i. e. if one of the 
bounding speeds 5^,5^ vanishes. 

A special case of HLLEM is obtained if we set = —S^ = |Sj„axl- djmamics and 

shallow water flow, the numerical viscosity of the resulting scheme coincides with the 
viscosity of the Rusanov/LLF (Local Lax Friedrichs) scheme for nonlinear waves and with 
the viscosity of the standard Roe scheme for shear and entropy waves. In the following, we 
refer to that scheme as LLFEM. 


3.2.4 HLLEMCC 

If in the definition of the HLLEM scheme, we replace the anti-diffusion coefficients 5^ 
by (1 — with 4> ^ [0,1], we can smoothly vary between HLLEM i<p = 0) and 
HLL i(j) = 1), i. e. between a complete and an incomplete Riemann solver. This is a 


technique which we used for our carbuncle cure of HLLEM (HLLEMCC) in [ ]2^ for the 
Euler equations and in 


for the shallow water equations. 

Since the goal of HLLEMCC is to prevent unphysical shear and entropy waves, its core is 
the computation of the residual in the Rankine-Hugoniot condition for these waves: 


t = fiQr) - fiQi) - u (^r - Qi) 


(19) 


with the Roe mean value for the normal velocity u. Now, we take as our basic indicator the 
residual relative to c, respectively its Euclidean norm 


e = ^llr 


12 5 


( 20 ) 


which vanishes for all shear and entropy waves. By introducing parameters a, /3 e (0,1), 
we can now define 


0(0, Fr^) = min{l, s 9 max{0, (1 —Fr“)}}^ 


( 21 ) 


for shallow water and 


0(0, MJ = min{l, (e Of max{0, (1 - M“)}} 


( 22 ) 


for the Euler equations, where F and are the directional Froude and Mach numbers 
perpendicular to the cell face. A major advantage of this approach is that the modification 
can be applied to shear and entropy waves separately. This we will use in Section 6.2 to 
somehow measure the influence of both wave types on the carbuncle. 

However, we still have to choose the parameters. Here, we adhere to the values as 
already published [[3,25,28,29]: For both shallow water and gas dynamics a = ^ = .33 
turned out to be a good choice. Furthermore, if not stated otherwise, we use e = 1/100 in 
the gas dynamics case and for the shallow water equations s = 10“^. 
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3.2.5 Shock fix 


While HLLEMCC modifies the numerical viscosity of HLLEM on the linearly degenerate 
waves, here we present a modification of the numerical viscosity at shocks. We will restrict to 
steady shocks in the Euler equations. Again we employ the residual in the Rankine-Hugoniot 
condition but this time on the nonlinear waves: 


= fiQr) - fiQi) - (d - c)(qr, -qi) = c[W 2 + W^ + 2 W 2 ) , (23) 

for shocks in the left wave, and for shocks in the right wave 

i^ = /(gr)-/(g!)-(“-c)(qr-gz) = -c(2^i + ^2 + ^3) ■ (24) 


Instead of using the residual in the Rankine-Hugoniot condition relative to the speed of 
sound, we employ the sum of the fluxes as our weighting factor. If we take the shock 
indicator 


9. 


1 - 


ll/(gr)-/(gz 


(25) 


V \\fiqr) + fiqi)\\2j 

with some small positive parameter we can modify the Einfeldt choice of wave speeds ( [T^ 
to 


S^=mm{eXui-Ci) + il-ejiu-c), u^-Ci, 0 }, 
Sjj = max{0s(“r + Cr) + (1 — 05 )(fi —c), u^ + c^, 0} 


resulting in a HLLEM-solver that does not enforce the Rankine-Hugoniot condition at steady 
shocks. In fact, the numerical viscosity is slightly increased as like in the Osher scheme, 
which we describe next. 


3.3 Osher-Solomon 


The Osher-Solomon scheme [44], a generalization of the Enquist-Osher scheme to 
systems of conservation laws, was designed with two major goals: prevent the sonic glitch 
and represent physically steady states as discrete steady states. For that purpose, they first 
had to abandon the Rankine-Hugoniot condition at shocks as a design principle. Instead, 
they use (generalized) Riemann-invariants. While this would lead to a slightly increased 
numerical viscosity at shocks, it also allows for stable steady state representations of steady 
shocks. Thus, the scheme is a good candidate for our purposes: measure the influence of 
the stability of the shock position on the carbuncle. The resulting numerical flux function 
reads as 

OosheriQu Qr) = ^ + /(^z)) “ X 


\Aiq)\ dq. 


(27) 


where T is a path in the state space which connects the left and right state via integral curves. 


Later on, Hemker and Spekreijse [22] presented an alternative choice for T, which they 


claim leads to better results. But, since in our tests we could not find any difference between 
the results of both versions, we just refer to it as Osher-scheme without differentiating for 


the integration path in (27). 
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4 Discussion of the test cases for the numericai 
investigation of the carbuncie 

In our discussion of the carbuncle and its sources, the choice of test cases plays an important 
role. They include classical examples like the double Mach reflection and the blunt body 
problem, which gave rise to the name of the carbuncle [|^|, as well as some -dimensional 
problems like the colliding flow problem, the steady shock, and the Quirk test. These tests 
have in common that they originate from one-dimensional problems which are artificially 
augmented with an additional space dimension. Since in these physically one-dimensional 
problems there is no inherent source for the carbuncle, we have to trigger it by adding 
some noise to the initial data. While most of these test cases have in common that in a 
physical sense we would not expect to see a carbuncle, the opposite is observed for the 


Elling test as drawn from the considerations in Section 2.3 


4.1 Blunt body problem 

In gas d 5 mamics, a popular test case for the carbuncle is the flow around a cylindrical 


obstacle |]9,12, M|d^d^d^||530 . In the shallow water case, a similar test would be the 
flow around a cylindrical bridge pier. We chose a pier with radius r = 1. As computational 
domain we employ a third of the annulus with inner radius r = 1 and outer radius R = 2. 
Since the interesting part of the flow is the inflow region, we restrict the domain in angular 
direction to The domain is discretized with 150 cells in the radial direction and 800 

cells in the angular direction. The initial flow is set to the inflow state everywhere. At the 
pier we employ wall boundary conditions, at the other boundaries first order extrapolation. 

Although, in principle, it is possible to do comparisons for blunt body flow also in the gas 
dynamics case p^ , here we restrict our investigation to shallow water. Since we want to 
start with the initial flow set to the inflow state, we could employ the gas dynamics version 
of the Osher scheme only for subsonic inflow. Thus, in this paper, we only show results for 
the shallow water case. For the inflow we choose Froude number Fr = 5. 


4.2 Double Mach reflection 

A famous test for the quality of a Riemann solver is the Double Mach reflection problem. It 


was suggested by Woodward and Colella [ 62 ] as a benchmark for Fuler codes. An analytical 


treatment is found in [39], O^, and references therein, while experimental results are 


presented in [ 19] and also in []^ pp. 152 and 168]. The problem consists of a shock front 
that hits a ramp which is inclined by 30 degrees. When the shock runs up the ramp, a 
self-similar shock structure with two triple points evolves. The situation is sketched out in 
Figure]^ To simplify the graphical representation, the coordinate system is aligned with 
the ramp-as done for the numerical tests. In the primary triple point, the incident shock i, 
the mach stem m, and the reflected shock r meet. In the double Mach configuration, the 
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Figure 2: Sketch of Double Mach Reflection with the carbuncle-t 5 ^e numerical artifact 
resulting in a kink in the leading Mach stem m. 


reflected shock breaks up forming a secondary triple point with the reflected shock r, a 
secondary (bowed) Mach stem m', and a secondary reflected shock r'. From both triple 
points, slip lines emanate. The reflected shock r' hits the primary slip line s causing a 
curled flow structure. 

While the main challenge for a high resolution scheme is to resolve the secondary slip 
line [3^ 62], the main challenge for first and second order schemes is to correctly 
represent the leading Mach stem. For some lower order schemes, the lower end of the Mach 
stem moves too fast, leading to a kink in the stem [48 ]. The situation is comparable to that 
for the blunt body problem. According to the wall boundary conditions, the ramp can be 
interpreted as the symmetry-line of a free-stream flow. Since the Mach stem is not perfectly 
aligned to the grid lines of a Cartesian grid, the same mechanisms come into place as for 
the blunt body problem. Thus, we expect Riemann solvers which produce a carbuncle for 
the blunt body to behave similar in that case. If the flow is split at the sjmimetry-line by 
reinterpreting it again as a wall (or ramp), we have to expect the lower part of the Mach 
stem to be kinked as sketched out in Figure 

For the numerical tests in this paper, we follow the guidelines in [30]. That means the 
boundary conditions at the upper boundary model a slightly smeared shock, and the vertical 
size of the computational domain is doubled compared to [62]. The only difference is that, 
instead of the vertical momentum, we choose the entropy for our plots. 


4.3 Colliding flow 

This test [ ]38] Section 7.7] resembles a simplified model for the starting process of the 
blunt body test when using the inflow state as initial data in the complete computational 
domain. This is best understood when considering the flow before the blunt body along 
the symmetry line. Since the flow is aligned with that symmetry line, and due to the switch 
of the sign of the flow velocity in wall boundary conditions, it behaves essentially like the 
left half of a colliding flow in Id. To turn it into a 2d-test, the flow is equipped with an 
additional space direction, in which everything is expected to be constant. In order to 
trigger the carbuncle, the initial state is superimposed with noise that is generated randomly 
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and has a small amplitude. 

For the gas dynamics test, in the initial state, density and pressure are set to p = 1, p = 1. 
The normal velocity is set to Uieft/right = ±20, the transverse velocity component to v = 0. 
To trigger the carbuncle, we superimpose artificial numerical noise of amplitude 10“^ onto 
the primitive variables instead of disturbing it in just one point as was done originally by 
LeVeque [ |3^ Section 7.7]. The computations are done on [0,60] x [0,30] discretized with 
60 X 30 grid cells. 

For shallow water tests, we set the initial height to h = 1 and the left and right velocities 
to u = ±30. The transverse velocity is zero. Onto this initial state, we superimpose artificial 
numerical noise of amplitude 10“^, but in this case we add the noise to the conserved 
variables. The computations are done on [—2.5,2.5] x [—2.5,2.5] discretized with 40 x 40 
grid cells. 

Since the problem is a simple 2d-extension of a one-dimensional problem, the results are 
presented in scatter-type plots: we slice the grid in x-direction along the cell faces and plot 
the entropy for all slices at once. 


4.4 Steady shock 

while the colliding flow test models the starting process of the blunt body flow, the steady 


shock test, introduced by Dumbser et al. [ 12 ], features a simplified model for the converged 
shock in the blunt body flow Following Dumbser et ah, we set in the upstream region p = 
1, u = 1. The upstream Mach number is set to M = 20, the transverse velocity component 
to V = 0. The shock is located directly on a cell face. To trigger the instability of the 
discrete shock profile, we add artificial numerical noise of amplitude 10“® to the primitive 
variables in the initial state. The computations are done on [0,100] x [0,40] discretized 
with 100 X 40 grid cells. 

For shallow water we set the water height and the Froude number at the inflow to h = 1 
and Fr = 30 respectively. At the outflow, we simply employ extrapolation boundary 
conditions. Again, we add artificial numerical noise, this time of amplitude 10“^, to 
the conserved variables in the initial state. The computations are done on [—2.5,2.5] x 
[—2.5,2.5] discretized with 100 x 40 grid cells. 

Again for the presentation of the results, we employ scatter-t 5 q)e plots as described for 
the colliding flow problem. 


4.5 Quirk test 


Quirk [48] introduced a test problem which is known as Quirk test. Contrary to the 
preceding example, it is not a one-dimensional Riemann problem, but consists of a shock 
running down a duct. The shock is caused by Dirichlet-tjqje boundary conditions on the 
left boundary with p = 5.26829268, u = 4.86111111, p = 29.88095238, while the flow 
field is initialized with p = l, u = v = 0, p = l/j- Originally, a disturbance of the middle 


grid line was used to trigger the instability [48]. Because the computations are done with 
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a Cartesian code, we instead use numerical noise in the same manner as for the steady 
shock and the colliding flow problem. The only difference lies in the amplitude of the 
perturbation, here 10“^. The computations are done on [0,1600] x [0,20] discretized with 
1600 X 20 grid cells. In this study, we only perform the test for the Euler equations. For a 
similar test in shallow water, we refer to [[^. 

Again we use scatter-type plots (as described above) to present the results. 

4.6 Elling test 

The Elling test is the experiment which was already described in Section The initial 
condition is a modified version of the steady shock test, cf. Section The region to the 
right of the shock remains unchanged. In the supersonic inflow region, only the middle 
x-slice is changed. Here the velocity is set to zero. This is done to model a vortex layer 
hitting the shock front. Again, we only perform the test for the Euler equations and refer 
to []^ for a similar test in shallow water. 


5 Influence of one-dimensional issues 

Although the carbuncle is a multi-dimensional issue, it is obvious that, especially on 
Cartesian grids, one-dimensional issues can drive multi-dimensional effects. 


5.1 Instability of shock position in first order schemes 

Two of the instabilities discussed in Section are purely one-dimensional. Both are 


instabilities of the shock: the instability of the shock position relative to the grid (Sec. 2.1) 


and the instability and ambiguity of the shock position within the grid cell (Sec. 2.2). 


Figure [^illustrates how the instability of the shock position can destroy a discrete shock 
profile. As Bultelle et al. [|^ point out, the shock position might jump by up to two grid 
cells. Figure [^ shows the worst case scenario: in one grid slice, it jumps two cells upstream, 
in the neighboring slice, it jumps two cells downstream. Thus, at a length of four grid faces. 


we created a new Riemann problem with all types of waves [28]. In the depicted situation. 


a strong flow downwards would be initiated. The instabilities and the ambiguity discussed 


in Section 2.2 can affect the discrete shock profile in a similar way. Since they are highly 


sensitive for perturbations, cross-flow might be induced within the grid slice containing the 
original shock itself. 

This situation could be avoided if the Id shock position was stable. However, as seen 
from the discussions in Sections 2.1 and 2.2[ this would mean to abandon the requirement 
of Riemann solver which exactly reproduces the Rankine-Hugoniot condition at a single 
shock. As mentioned above, the Osher solver replaces this by the requirement of yielding 
steady discrete solutions for any steady discontinuity and employs Riemann invariants 
over all nonlinear waves. Thus, the Osher scheme seems to be a good candidate to avoid 
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original shock location 



Figure 3: Unstable two-dimensional shock profile 
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Figure 4: Scatter-t 5 ^e plots of entropy for steady shock with different numerical fluxes 
at t = 1000. 
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Figure 5: Scatter-type plots of entropy for Quirk test with different numerical fluxes at t 
125. 


the instability of the Id shock position, much the same as the HLLEM-scheme with the 
modification at steady shocks (Sec. 3.2.5). Another solver which also disregards the 
Rankine-Hugoniot condition is LLFEM. The difference to Osher and the steady shock fix 
is that the numerical viscosity on the nonlinear waves is much higher. As the results in 
Eigure show, the LLEEM cannot prevent the carbuncle while the other two can. Erom 
that, we conclude that it is important to add the correct amount of viscosity in order to 
stabilize the shock position. 

Figure [^reveals that if we guarantee steady discrete representations of steady states this 
is not sufficient to guarantee a proper representation of unsteady flows. The Quirk test is 
employed and shows that the Osher scheme still might produce carbuncle-like structures, 
although much weaker than with, e. g., HLLEM. 

Our setting for the blunt body problem results in a highly transient starting phase which 
eventually passes into a steady state. This raises the question which property of the Osher 
scheme would be dominant: the tendency to produce a carbuncle in a transient flow or the 
guarantee for steady discrete representations of steady flows. Thus, in Figure]^ we show 
the shallow water flow around a cylindrical pier at different times. At time t = 1.5, it is 
easy to see that the scheme produces a slight carbuncle-type structure. When the discrete 
flow field eventually becomes steady, obviously that carbuncle is somehow smoothed out. 
The solution is steady but not physical, confirming the statement by Robinet et al. [49 ] that 
the carbuncle can lead to unphysical steady states. 


5.2 Influence of the order of the scheme 

Since it is often reported that increasing the order of the scheme applied in the computation 
reduces the carbuncle, here we investigate the relation between the order of the scheme and 
the Id-stability of the discretized shock. At this point, we should stress out that throughout 
this paper the term order refers to the design order of the scheme, which is only achieved 
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Osher at t=0.75 


Osher at t=1.5 


Osher at t=2.25 


Osher at t=3.0 






Figure 6: Flow around cylindrical pier in shallow water with 1st order Osher scheme at 
different times. Energy shown. 


in smooth parts of the flow field, and not to the actual order of the scheme which would 


automatically drop in the vicinity of shocks. As was pointed out by Roe and Zaide [63 -65 ], 
a major role is played by the fact that it is impossible, at least for gas dynamics and shallow 
water flows, to split a shock satisfying the Rankine-Hugoniot condition into two consecutive 
shocks which both would satisfy the Rankine-Hugoniot condition. This situation improves 
for higher order, at least when geometric reconstruction is applied. For the computation of 
the inter-cell fluxes, two states are available, one at the left cell boundary and one at the 
right cell boundary. Thus, for a steady shock it would be sufficient to satisfy 


Kqi) = fiqj = = fiqr), 

which can easily be achieved by 


(28) 


q^ = qi, qt = qr- ( 29 ) 

This is (for a scalar situation) sketched out in Figure]^ For the sake of simplicity, in the 
following, we restrict our considerations to the scalar case. It is easy to see how the results 
can be transferred back to the systems case. 

Figure]^ shows how the situation improves with increasing order of the scheme. If, for 
instance, we employ a polynomial reconstruction in the cell where the shock is located, for 
high orders, condition ( |2^ can be ensured for almost all shock positions without sacrihcing 
monotonicity: 


18 












1st orde 

r 

^r 

ti order 







2nd ord 
very hig 

1 


qr 











Qi 



q q~^ 

“m 



-► 


X 

Figure 7: Reconstruction of shock profile (dashed line) on finite volume grid with different 
orders 


Theorem 1. Given the situation depicted in Figure]^ let the cell with the shock be ■^i+ 1 / 2 ] 

and Ax = Xj +^/2 ~ ■^i-i/ 2 - Furthermore let Aq = q,. — qi 7 ^ 0. Let the shock be located 
at Xj_i /2 + 0 Ax with 9 e [0,1]. 

1. For any 0 e [0, a polynomial reconstruction with degree less or equal two can 
be found such that condition ( [ 2 ^ is satisfied. 

2. In the cell [Xj_i /25 ^ 1 + 1 / 2 ! reconstruction can be made monotone for polynomial 

degree <nif 9 ^]. 

Proof. Statement!^ is obvious and well known. It was already used by van Leer in his work 
on higher order methods []60j]. 

For the proof of statement we can assume, without restriction, 

[Xi_i/ 2 ,Xi+i/ 2 ] = [ 0 , 1 ], qi = 0, qr = l, and 0 > 1/2 . 

All other cases can be derived from this by symmetry, scaling, and translation. 

It is easy to see that with above settings 

1 

q(x) dx = 1 — 0 . (30) 

Jo 

On the other hand, we know for any monomial 

1 

x"dx =-. (31) 

Jo ^ + 1 

Furthermore all monomials with degree greater of equal one are monotonously increasing 
in [0,1] and attain the values 0 at x = 0 and 1 at x = 1. Obviously, the same is true for 
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all weighted means of such poljmomials as long as the weights are non-negative. Thus, 
for any 6 e [|, we can find a polynomial p„ of degree < n, monotonously increasing 
in [0,1], with 


Pn(0) = 0, p„(l) = l, 


and 


Pnix) dx = 1 — 9 . 


(32) 


As already mentioned, the general case follows from this by symmetry, scaling, and transla¬ 
tion. □ 

Note that the linear reconstruction of case in Theorem is not achieved by standard 
second order schemes. Due to their restricted stencil, they cannot distinguish between the 
situation of Theoremcase(slope Aq/ Ax in the middle cell), and a linear state with 
slope Aq/{2Ax). But in order to achieve second order, they have to reconstruct linear 
states exactly. Slope Aq/Ax in the middle cell can be achieved, e. g., by the limiter cpswebyj 
the upper bound of the Sweby region, as described in [26, Section 2.4.2], leading to a 
first order scheme. For other 9 e [—1,1] it would still satisfy one of the identities ( |2^ . 
While Minmod cannot satisfy any of the identities ([2^ for 9 e (—1,1), Superbee shares 


the behavior of Pswehy for 1^1 ^ 1/3 and the MC-Limiter for |0| > 1/2. 

In this context, it is also worth to note that some third order schemes, although formally 
based on linear reconstructions, e. g. [[^[^ in the final analysis still employ parabolic recon¬ 
structions. For each cell, they compute two linear reconstructions in the manner described 


in [26 Section 2.4.2], one ensuring third order for left-going waves and another one 


ensuring third order for right-going waves. While is taken from the first reconstruction, 
q/^ is taken from the latter. Although it is, in general, not possible to reinterpret this as 
a single linear reconstruction, it is always possible to reinterpret it as a parabolic recon¬ 
struction. But due to their restricted stencil, they do, in general, not satisfy the identities 
in equation ( |2^ . Strangely enough, none of the authors of such limiters considers the 
resulting parabolic reconstruction, and, therefore, none of them checks if this parabolic 
reconstruction is indeed monotone. 

Since many modern schemes, like ENO/WENO, do not explicitly enforce monotone 
reconstructions in the cells, for these schemes, we can expect identities ( |2^ to be satisfied 
in most cases, and, if not so, at least approximated with very small error. Thus, for 
schemes of order greater or equal three, we we might expect physically steady shocks to be 
represented as discrete steady shock, independent of the Riemann solver and the position 
of the shock, and hope for a similar behavior for moving shocks. 


6 Influence of two- or three-dimensional issues 

For years, the research on the carbuncle concentrated on two-dimensional issues, although 
the Osher scheme was already in widespread use, and the carbuncle was known to be 
rarely found in very high order schemes or on unstructured grids. Here, we started with 
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Figure 8: Steady shock problem in shallow water at times t = 0.2, 0.4, 0.6, 0.8, 1.0 with 
Roe scheme, energy shown. Upper row: first order; second row: second order 
with Superpower limiter. 



Figure 9: Colliding flow problem in shallow water at times t = 2/3, 1, 4/3, 5/3, 2 with 
Roe scheme, energy shown. Upper row: first order; second row: second order 
with Superpower limiter. 
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Figure 10: Leading shock structure of double Mach reflection in gas d 5 niamics at t = 0.2 
with different numerical fluxes, entropy shown. Upper row: first order, lower 
row: second order. 
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Figure 11: Flow around cylindrical pier in shallow water at Froude number Fr = 5 with 
Roe and Osher, first order and second order with Albada 3. 
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Figure 12: Scatter-type plots of entropy for Quirk test with different numerical fluxes and 
second order at t = 125. 
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Figure 13: Transverse velocity at different times in gas dynamics steady shock problem for 
different solvers. 


one-dimensional issues to better understand how they, when everything is generalized to 
two or three space-dimensions, interact with two- and three-dimensional issues. Hence, in 
this section, we want to concentrate on the interplay of Id and 2d (or 3d) issues. 


6.1 Numerical shear viscosity 


At this point, we have to refer the reader back to Figure It is obvious that the vertical 
flow induced by the instability of the shock position in turn induces a strong shear flow 
not only in the newly created Riemann problem, but even more along the original shock 
profile. In summary, some kind of turbulence at the original shock position is created 
which is superimposed onto the original flow. Due to its construction, the HLLEMCC solver 
distinguishes between shear waves which are superimposed onto nonlinear waves and 
shear waves which are not. On the first, it behaves like HLLE, thus damping the turbulence, 
on the latter, it behaves like HLLEM, here allowing, e. g., for well resolved boundary layers. 
In our previous works [j^25 28 2^, we could show that HLLEMCC prevents the carbuncle 


while allowing for good resolution of physical shear waves. 

In Eigure[T^ we demonstrate how the above described mechanism drives the carbuncle. 
We show the transverse velocity at different times for HLLEM, HLLE, and the Osher scheme. 
By stabilizing the Id shock position, the Osher scheme keeps the transverse velocity at 
about the magnitude of the artificial numerical noise introduced in the initial state. Thus, 
no carbuncle arises. The HLLE scheme, by its excessive shear viscosity, keeps the transverse 
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Figure 14: Quirk test at t = 125 with 1st order standard HLLEMCC and different versions 
of 2nd order HLLEMCC scheme. 
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Eigure 15: Steady shock in gas dynamics at t = 1000 with 1st order standard HLLEMCC 
and different versions of 2nd order HLLEMCC scheme. 


velocity below some threshold and thus also avoids the carbuncle. For HLLEM, there is 
no mechanism to damp the turbulence along the original shock. Over time, a carbuncle 
evolves. 


6.1.1 Shear viscosity and order of the scheme 


As we have seen in Section [5^ higher order has a stabilizing effect on the Id shock position. 
This raises the question if for second order schemes, the carbuncle correction in HLLEMCC 
might be relaxed. The answer is not obvious since raising the order also lowers the viscosity 


on the shear waves. In Figure 14 we give a comparison of first order HLLEMCC and 


several implementations with second order. The difference between the versions is in 


the choice of the parameters. As mentioned in Section 10 for the standard HLLEMCC, 


the parameter e in equations pT] ) and ( [22] ) is chosen as £ = 0.01. Here we also show 
results for e = 0.005, e = 0.00125, and for s = 10“^. For the latter, the results are almost 
indistinguishable from the pure second order HLLEM. All in all, the results suggest to leave 
the parameters unchanged and stay with the same parameters as in the standard version for 
first order. The gain in stability of the Id shock position and the loss in shear viscosity are 
just in balance. Eor the steady shock test in Eigure [Tsj the situation slightly improves. But 


it is still recommended to use HLLEMCC with the set of parameters given in Section 3.2.4 


6.1.2 Shear viscosity and the resolution of physical carbuncles 

In Section [23 we discussed the findings of Elling [|T^ on physical carbuncles. This lead us 


to the Elling test case as described in Section 4.6 which allows us to test the numerical flux 


functions for their ability to resolve these physical carbuncles correctly. As we can see from 
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Figure 16: Elling test for gas djuiamics at time t = 100; comparison of different solvers; 
entropy shown. 


Figure the numerical results with Osher, HLLEM, and HLLEMCC are rather similar even 
with the first order scheme. The HLLE scheme, which suppresses the carbuncle by a severe 
amount of numerical shear viscosity, tries to prevent even the physical carbuncle, not only 
in the first order, but also in the second order computation. Since the test case is a model 
for shock boundary-layer and other shock vortex interactions, we note that schemes which 
rely on the HLLE flux, even when it is only applied locally in the vicinity of strong shocks, 
might destroy some physical features of the flow. Thus, codes which are based on a switch 
between complete and incomplete Riemann solvers depending on the distance to the next 
strong shock should be carefully tested with the Elling test before applying them to more 
complex flow problems. If the physical carbuncle is not properly reproduced, the switching 
mechanism has to be reworked. 


6.2 Influence of viscosity on entropy waves 


Some authors consider the carbuncle a result of the treatment of mass transport and entropy 
waves [40] within the Riemann solver. The hunt for entropy consistent Riemann solvers, 
e. g., is at least partially driven by that idea. And indeed, for some schemes this causes 
problems, e. g. for Elux Vector Splitting (EVS) schemes, which may loose positivity by 


exactly resolving entropy waves [ 18 ]. But there is no strict proof, not even the proof by 


Liou and Steffen in [40], for a connection between the resolution of entropy waves and the 
carbuncle. 

From our considerations in Section]^ we know that the instability of the Id shock position 
causes a new Riemann problem perpendicular to the original shock, which includes all 
types of waves (for gas dynamics also entropy waves). Thus, although the carbuncle occurs 
also in shallow water, where there are no entropy waves, we can conclude that there is 
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Figure 17: Steady shock in gas d 5 mamics; comparison of HLLEM without carbuncle cor¬ 
rection, correction on entropy waves, and full HLLEMCC at t = 1000; entropy 
shown. 
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Figure 18: Colliding flow in gas d 5 mamics; comparison of HLLEM without carbuncle cor¬ 
rection, correction on entropy waves, and full HLLEMCC at t = 20; entropy 
shown. 


a connection between carbuncle and entropy waves. The question we have to answer is: 
Which type of waves has the stronger impact on the stability of discrete shock profiles: 
shear waves or entropy waves? 

A good means to answer that question is the HLLEMCC solver, which allows to apply the 
carbuncle correction to both types of linearly degenerate waves separately. In Figures [Tt] 
and[^ we present results for the steady shock and the colliding flow problem. We compare 
pure HLLEM with three versions of HLLEMCC: correction only applied to entropy waves, 
correction only applied to shear waves, full HLLEMCC. Erom the numerical results we 
easily conclude that the resolution of entropy waves contributes to the carbuncle, but the 
contribution is small compared to the contribution by the shear waves. 


7 Conclusions and directions for further research 

In this paper, we investigated the origin of the carbuncle phenomenon. Guided by the 
theoretical results reviewed in Section we found a set of numerical test cases which 
helped us to sort out the different issues involved in the carbuncle. We observed that the 
Osher scheme by its special construction of the numerical viscosity on shocks suppresses 
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the carbuncle to a certain extent. For steady grid-aligned shocks, we could remodel this 
viscosity in the HLLEM scheme, confirming that it is in fact the numerical viscosity on the 
shock which is responsible for the stabilizing effect. We also showed how increasing the 
order offers an alternative way to stabilize the shock position: introducing more degrees of 
freedom for the scheme allows to remodel the Rankine-Hugoniot condition at a captured 
shock. This stabilizing mechanism even works when the mechanism of the Osher scheme is 
not sufficient an 5 anore: in the case of not perfectly-grid aligned shocks like the bow shock 
in the blunt body problem. 

Since the instability of the Id shock position creates vorticity along the shock, we also 
considered the influence of the numerical viscosity on entropy and shear waves. We found 
that the influence of the shear viscosity is much higher than that of the viscosity on entropy 
waves. We also found (cf. Section [2^ that Riemann solvers like HLLEMCC, which try to 
reduce the carbuncle from within the Riemann solver and without too much loss of the 
resolution of shear layers, should not be altered when used in a higher order scheme. The 
gain in stability of the shock position is compensated by the reduction of the shear viscosity. 
By employing the Elling test we could show that incomplete Riemann solvers like HLLE 
not only prevent non-physical carbuncles but also the physically induced breakdown of the 
shock profile when it is hit by a vortex layer. 

What is still lacking, is a deeper understanding of the amount of numerical viscosity 
on shocks in order to stabilize the shock position already for low order schemes. But 
this would be desirable when one wants to combine the stabilizing mechanisms of Osher 
and HLLEMCC, in which case the shear viscosity of HLLEMCC could be further reduced. 
Furthermore, one would hope to be able to extend the theoretical results reviewed in 
Sectionj^to the case of higher order schemes. As we have seen, there is a significant impact 
of the order of the scheme on the stability of discrete shock profiles. 
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